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ABSTRACT 




A technique for reducing the number of iterations necessary 
for solving linear programs using the primal-dual algorithm is 
presented. It appears that the new method will also decrease the 
number of iterations over any other simplex algorithm. A FORTRAN 
program incorporating the technique, as well as some comparative 

X 

computational results are given. 
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SECTION I 



INTRODUCTION 

The primal-dual algorithm, as developed by Ford and Fulkerson, 
theoretically permitted faster convergence in solving linear programs 
over such methods as the two-phase or Charnes 1 M-method . Instead of 
first driving the artificial variables to zero to obtain a feasible 
solution and then working to optimality in a second phase, the primal- 
dual works towards optimality and feasibility in the primal and dual 
problems at the same time. This is accomplished during what would 
be the first phase for the other procedures. However, in actual 
computational experience, the primal-dual algorithm has not achieved 
faster convergence over the other simplex procedures. 

This paper presents^* a modification of the primal-dual algorithm 
that does accomplish faster convergence than the usual primal-dual 
algorithm and perhaps also the other simplex methods. 



Professor Harold Greenberg provided the idea for the new 
technique. 
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SECTION II 



THE STANDARD PRIMAL- DUAL ALGORITHM 



Using matrix notation, the linear programming problem is: 

Minimize z = z + CX 
o 

Given AX = H 

x^^O j = 1 , . . . ,n 
where z is a constant and 



C is a 1 x n vector with components c_. ^ 0, j = l,...,n 

H is an m x 1 vector with components h^ ^ 0, i = 1 m 

A is an m x n matrix with elements a. . i = l,...,m 

i > J > 

j = 1 » • • • > n 

X is an n x 1 vector with components x., j = l,...,n 

2 

The standard primal-dual, as it appears in Dantzig, consists 



in augmenting the constraint equations with artificial variables as, 



AX 4- Y = H 

^ 0 i = 1 m 



( 1 ) 



where Y is an m x 1 vector with components consisting of the arti- 
ficial variables y^; and forming 



w = w + DX 
o 



where 



w 

o 



m 

= E 
i=l 



h. 

l 



and D is a 1 x n vector with components d., j = l,...,n where 



George B. Dantzig, Linear Programming and Extensions 
(Princeton: Princeton University Press, 1963), pp. 247-53. 
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d. 

J 



m 




It is noted that both the c. and the h. must be positive 

J i 

prior to starting the computations. Negative tu may be made 
positive by multiplying the associated constraint equations by 
a minus one. The method for correcting negative c^ is discussed 
in the "Description of the Computer Program" section of this 
paper . 

The problem is to find x_. ^ 0, w = 0, that minimizes z 
subject to (1). At any stage of the iterative procedure the 
equations have the form: 

X B + AjX = 

Y + A 2 X - H 2 (2) 

DX = w - w 

o 

CX = z - z 

o 

where X is a vector representing the current basic variables 

D 

and X is a vector representing the current non-basic variables. 

The numerical subscripts 1 and 2 represent partitioned matrices 
and vectors. Matrices and vectors with bars are the corresponding 
original matrices and vectors after the usual pivot operations 
within the primal and dual problems. When an artificial variable 
becomes non-basic, it is dropped from further consideration. 

This results in the matrix and the vector. Each iteration 
in the standard primal-dual algorithm consists of the following: 
a) Finding the most negative d^ value for those indices 
j having "c. =0. 
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b) If there are no indices with this property, the 



multiplier k = min [c ./ - d.] is used to obtain at 



d .<0 J 
J 



least one c. equal to zero by 



c . = c . + k(d . ) 
J J J 




c . = 0. 
J 



Then x is selected to be basic, 
r 



d) Then min (h./a. ) is found for a. > 0. This will 

l i,r i,r 



select a pivot row in either the X or Y sets of 

D 

equations . 

e) The simplex pivot method is then used. If the pivot 

row is in the Y set of equations the artificial variable 
in that row is dropped. Thus, in any case, one of the 
forms in (2) is achieved. 

When w = 0, the basic solution is optimal and the problem 
is solved. 

The actual computations should include rules to prevent 
cycling, such as the perturbation technique or the use of 
lexicographic ordering. 
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SECTION III 



DESCRIPTION OF THE NEW TECHNIQUE 

The primal-dual algorithm is modified in a two-step procedure. 

In the first step, a pivot row in the Y set of equations is selected. 
By only considering elimination of artificial variables, after m 
or less iterations, w equals zero. At this point the equations 
have the form: 

X + AX = H 
B 

z + CX = z 
o 

where all c^ are positive. If all the tu are positive, the solution 
is optimal. We then stop. However, in most cases, this stage of 
the problem will contain one or more negative h^. In this case 
the subscript "it" is assigned to the row with the most negative 
tu. The following procedure is then used as the second step of the 
modification: 

a) The equation containing h^ t is multiplied through by minus 

one. This produces h. >0. 

it 

b) The "it” equation is then added to all equations with 
negative tu. All basic variables except the "it" equation 
now become feasible. 

c) To maintain the canonical form an artificial variable is 
added to the M it M equation. A new infeasibility equation 
with w = h. is formed from the "it" equation. 

The process now reverts to the primal-dual algorithm to achieve 
w equal to zero. Since the canonical form has been maintained and 
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all c have been kept positive throughout the process, 
w equals zero, the solution is optimal. 



once 
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SECTION IV 



COMPUTATIONAL RESULTS 

Several problems were run on an IBM 360/67 using both the 
standard primal-dual and the primal-dual with the new technique. 
The problems used were test problems from the SHARE Linear 
Programming Project and two highly degenerate assignment type 
problems. Results are given in Table I with the two assignment 
problems listed first. 



TABLE I 

TABLE OF COMPARATIVE RESULTS 



PROBLEM 


SIZE 
(m x 


n) 


SIMPLEX 

ITERATIONS'^ 


PRIMAL-DUAL WITH 
THE NEW TECHNIQUE 
Iterations Seconds 


STANDARD 
PRIMAL-DUAL 
Iterations Seconds 


Assignment 


74 


X 


150 




89 


32 


83 


29 


Assignment 


34 


X 


376 




74 


12 


97 


28 


SHARE 


1A 


33 


X 


64 


70 


39 


7 


39 


7 


SHARE 


ID 


27 


X 


45 


61 


38 




38 




SHARE 


IE 


31 


X 


106 


66 


96 


30 


136 


42 


SHARE 


IF 


66 


X 


135 


184 


125 


86 


173 


114 


SHARE 


2A 


30 


X 


103 


114 


79 


29 


119 


43 



In the two assignment problems the primal-dual with the new technique 
resulted in optima different from that of the standard primal-dual since 



Philip Wolfe and Leola Cutler, Experiments in Linear Programming,” 
Recent Advances in Mathematical Programming (New York: McGraw-Hill 

Book Co., Inc., 1963), p. 198. 
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these problems have non-unique solutions. The two methods resulted 
in the same optima in all other problems. 

No accurate timing data is available. The timing data in the 
above table though does give an indication of the reduction in 
execution times possible with the new technique. However, the 
modified method does perform less work than the usual method 
because in the first step ratios and comparisons are only taken 
for the artificial rows. The restart of the problem in the second 
step requires only an insignificant amount of time. Thus, in the 
table, the number of iterations represents a measure of improvement 
of the new method. The new technique shows significant reductions 
in the number of iterations from the two-phase standard simplex 
method in four out of the five problems where a comparison is 
possible. When comparing the new technique with the standard 
primal-dual, a similar reduction is evident in the larger problems. 
It therefore appears that the new method presents a significant 
improvement when the number of variables is large. 
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SECTION V 



DESCRIPTION OF THE COMPUTER PROGRAM 



The accompanying program (see Appendix) is written in the 
FORTRAN IV language, and was used to run all test problems with 
the exception of the two assignment type problems (since these 
only deal with zeros or ones in the A matrix, considerable 
reduction in the core memory requirement could be made with a 
slight modification of the program). The program will handle up 
to 69 constraint equations and up to 149 variables, and requires 
about 160,000 bits of core storage. Columns may be designated 
by alphanumeric names, but rows must be consecutive integer 
values not to exceed the value 69. Only the non-zero elements 
of A, C, and H need be read in. All h^ must be positive, but 
negative c^ are allowed. 

Names used in this program conform to those already outlined 
in this paper. Those not previously discussed are defined as: 



FF is a multiplicative constant used to insure that all a. 



i, J 



are integer. 

EP is a constant used in the round-off error control. 

DI is the cumulative multiplication of the pivot elements. 
CE is a vector identical to the initial C, and is used to 
calculate z. 

AX is the pivot column in the Revised Simplex tableau. 
CNAME is a vector of column names. 

DUAL is the vector of dual variables. 
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DB is a vector of simplex multipliers used to calculate the d ^ . 

R is a vector used to record whether each variable is basic 
or artificial. 

E is a computed scalar used in the updating and round-off 
control . 

P is a vector of multipliers used in the actual pivot operation. 

ITR indicates the number of iterations performed. 

INDEX controls the use of the new technique (INDEX equals zero 
yields the standard primal-dual, and INDEX equals one 
incorporates the new technique) . 

IB is a vector of variable designators (a negative number 

indicates an artificial variable, and a positive number 
indicates a basic variable). 

All other names are defined in the program. 

Once the data has been read in, the c_. are tested for negative 
values. If all c^ are positive, the primal-dual commences. However, 
if there is one or more negative c_. the following procedure is fol- 
lowed : 



a) RMIN is set equal to the minimum c.. 

J 

b) n is set equal to n + 1 and m is set equal to m + 1. 

c) The coefficients of all variables in the added constraint 
equation indicated in b) above are set equal to FF including 
x^. (Variable x^ does not appear in any other constraint 
equation.) 

d) Then c^ is set equal to minus RMIN. 
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e) The right hand side value, h^, is entered as: 

h = k(FF) 
m 

where k is a number larger than the sum of the 
in the optimal tableau. 

The problem is now modified so it may be handled by use of the 
primal-dual algorithm. 

The remainder of the program is a FORTRAN representation of 
the steps outlined in the sections "The Standard Primal-Dual 
Algorithm" and "Description of the New Technique," However, for 
computer calculations, an essential addition to the procedures is 
some control over round-off error. 

Round-off error develops when two large numbers, say y and z, 
(differing only slightly from each other), are subtracted from 
each other such as : 

x = y - z 

In this program the absolute value of x divided by y is compared 
to EP. If it is less than EP, x is set equal to zero. As a 
further control, the values of DB, B, and H are updated each 
iteration. 
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1 

2 



PROGRAM TO SOLVE PRIMAL-DUAL LINEAR 
PROGRAMMING PROBLEMS 

IMPLICIT REAL*8 (A-H,0-Z) ,, . . , 

COMMON A (70, 150) ,B(70,70) ,C( 150) ,0( 150) ,H(70) , AX( 70 
COMMON DB (70 ,OUAL <70),R(150),CE(150),E,W,Z,DI,FF,EP 
COMMON CN AM E ( 150) 

COMMON 18(70) »M,N,JS, I R, INDEX, ITR,M1,N1 
F E A L * 6 NINE/' 9999 •/ 

FF=10000. 

EP=0 . 000C01 
0 1 = I . 

M=0 

DO 1 1=1,70 
H ( I ) =0 . 

DO 1 J=1 , l 50 
A ( I, J)=0. 

DO 2 1=1,150 
CNAME( I )=C. 

CE ( I ) = C • 

C( I ) =0. 



READING IN COEFICIENTS OF CONSTRAINT EQUATIONS 

DO 5 N= 1,1000 

PEAD (5,3) CB, J, RB 

IF(CB.EQ.NINE) GO TO 16 

3 FORMAT ( 6X , A6, 16, F10.4) 

DO 4 1=1,150 

IF(CNAME( I) .EQ.CB) GO TO 30 
IF(CNAME( I) .EO.O. ) GO TO 30 
GO TO 4 

30 A ( J , I ) =RB*FF 
CNAME ( I ) =CB 

WRITE (6,102) I, J, A ( J , I ) 

102 FORMAT (IX, 2110, F18.9) 

IF(J.LE.M) GO TO 5 
M=J 

GO TO 5 

4 CONTINUE 

5 CONTINUE 

READING IN RIGHT HAND SIDE VALUES 

16 DO 6 N=1 , ICO 
READ (5,14) I, RB 

14 FORMAT (118, F10.4) 

IF ( I ,GT ,999 ) GO TO 17 
H( I )=RB*FF 

WRITE (6,14) I, H ( I ) 

6 CONTINUE 

READING IN COSTS 

17 DO 9 J= 1 , 200 

PEAD (5,7) CC, II, RB 

7 FORMAT ( 6X » A6, 16, F10.4) 

IF(CC.EQ.NINE) GO TO 18 

DO 11 K=1 ,150 

IF(CC.EQ.CNAME(K) ) GO TO 12 
GO TO 11 
12 C(K)=RB*FF 
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C E ( K ) = R8 

WRITE (6,13) K, II, C ( K 1 
13 FORMAT ( 6X , 216, F15.5) 

GO TO 9 
11 CONTINUE 
9 CONTINUE 
18 CONTINUE 

DO 20 1=1,150 

IF(CNAME( I) .GT.C. ) GO TO 20 
N= I- 1 
GO TO 21 

20 CONTINUE 

21 WRITE (6,101) M, N 

101' FORMAT (10X, • NUMBER OF ROWS=', 15, 15X, 

1 'NUMBER OF COLUMNS=', 15) 

WRITE (6,103) (CNAME(I), 1=1,175) 

103 FORMAT < IX, 20A6) 

WRITE (6,104) (C(J), J= 1 , N ) 

104 FORMAT (IX, 10E12.4) 

I NDE X= 1 

N1 = N+ 1 
M 1 = M + 1 
CALL LIP 
8 STOP 
END 



21 



no ° non noon on 



SUBROUTINE LIP 



IMPLICIT RFAL*8 (A-H,0-Z) 

COMMON AI 7C, 15G) ,B(70,70) ,C < 150) ,D( 150) ,HI 7C > , AXI 70) 
COMMON l)B ( 7C ) t DUAL ( 70 ),R(150),CFI150>,E,W,Z,0I,FF,EP 
COMMON CNAM F ( 150) 

COMMON IB(70) , M , N , JS » IR , INDEX , I TR , Ml f N1 
E = 0. 5/D I 

MODIFYING THE PROBLEM IF THE OBJECTIVE 
FUNCTION CONTAINS NEGATIVE COSTS 



50 

52 



53 



RM I N=0 • 

DO 50 1=1, N 

I F ( C ( I ) . GE . E ) GO TO 50 
I F ( C ( I ) .GE.RMIN) GO TO 50 
RM I N=C ( I ) 

CONTINUE 

IFIRMIN.EQ.O.) GO TO 60 
DO 52 1=1, N 
C ( I ) =C ( I ) -RM I N 
N=N+ 1 
M=M+ 1 

HIM) = 1000 ♦ * FF 
C ( N) =-RMI N 
DO 53 J=1,N 
AIM, J)=FF 



INITIALIZING THE TABLEAU 

60 DO 3 J= 1 , M 
DO 2 1=1, M 

2 B I I , J ) = 0. 

3 B I J , J ) =1 . 

DO 4 1=1, M 
DB I I ) =- 1 . 

DUAL I I ) =0 . 

4 IBI I ) =- I 
DO 5 J= 1 , N 
R I J) =C. 

0 1 J ) =0 • 

DO 5 1 = 1, M 

D< J ) =0 I J ) -A ( I ,J) 

5 CONTINUE 
W=0. 

DO 6 1=1, M 

6 W=W + H I I ) 

I TR=0 



OBTAINING A COST EQUAL TO ZERO AND A PIVOT COLUMN 



7 

3 

9 

10 



IFIDABSI W)-E) 100,100,8 

DO 9 J= 1 , N 

IFIDI J) ) 1 0 , 9 , 9 

CONTINUE 

GO TO 101 

JS=0 



DO 11 J=1 , N 

IFIRI J) .EQ. 1. ) GO TO 11 
IFIDI J) .GE.-E) GO TO 11 
IF(JS.EO.O) GO TO 13 
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12 IF(CX+C(J)/D(J) )11,11,13 

13 CX=-C( J) /D( J) 

JS=J 

11 CONTINUE 
47 DO 14 J = 1 » N 

14 C ( J) =C ( J ) +CX*0( J ) 

C( JS)=0. 

DO 301 J = 1 » N 
IF(C(J).LT.E> GO TO 301 
GO TO 302 

301 CONTINUE 
GO TO 100 

302 Z=0. 

• DO 303 1=1, M 
J= I B ( I ) 

IF(J.LE.C) GO TO 303 
Z=Z+CE( J) *H( I » 

303 CONTINUE 

DO 40 1=1, M 

40 DUAL < I )=DUAL( I )-CX*DB ( I ) 

15 IF(W-E) 100, 100,41 

41 DO 16 J = 1,N 

IF (C(J 1-0.00001) 17, 17, 16 

17 C(J)=0. 

18 I F ( D ( J)-D( JS) ) 19, 16, 16 

19 JS=J 

16 CONTINUE 

WRITE (6,213) JS, DUS) 

213 FORMAT ( 1H, *JS=', 14, 3X, • D= • , F12.4) 
I F ( D ( JS) ) 20 , 7 , 7 

20 CALL ELEM 

IF ( D I - 1 . )528,26,26 
26 CALL PIVOT 
GO TO 15 

PRINTING OF RESULTS 



100 WRITE (6,204) 

204 FORMAT ( 1H , 'RESULT ’ ) 

GO TO 425 

101 WRITE (6,205) 

205 FORMAT { 1H , 'INFEASIBLE') 

STOP 

425 DO 426 1 = 1, M 

426 WRITE (6,526) IB(I), H(I), OUAL(I) 

526 FORMAT (1H , 'BASIS VARIABLE NUMBER= ' , 14, 

1 • VALUE OF BASIS VARIABLE^, F12.4, 5X, 

1 • DUAL= ' , F12.4) 

I F ( RM IN. E0» 0. ) GO TO 427 
ZZ=0. 

DO 429 K= 1 , M 
J= I B ( K ) 

IF(J.LE.O) GO TO 429 
ZZ = Z Z + CE { J ) *H( K ) 

CONTINUE 
GO TO 428 

11 = 1 

WRITE (6,527) W, ZZ 

FORMAT (1H , • W= ' » F12.8, 5X, *Z=', F12.5) 
IF(INDEX.EO.l) GO TO 1 
RETURN 



429 

427 

428 

527 

528 



5X, 
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THE NEW TECHNIQUE 



1 I T=0 

DO 30 1=1, M 

I F ( H ( I) .GT.-E) GO TO 30 
I F ( IT.EO.O) GO TO 31 
IF(H( IKGE.HI IT) ) GO TO 30 
31 I T = I 
30 CONTINUE 

IF( IT.EO.O) GO TO 528 
i e ( i t i =- 1 t 
HI I T ) =-H ( IT) 

DO 34 K=1 , M 
D8(K)=B( I T , K ) 

34 B ( IT » K ) =-B ( I T , K ) 

DO 33 1=1, M 

IFIHI I ) .GT.-E ) GO TO 33 
HI I ) =H I I ) +H I IT) 

DO 35 K= 1 , M 

35 B I l » K ) =B I I , K ) +B I IT,K) 

33 CONTINUE 

W=H ( IT) 

DO 36 J=L,N 
AX I I T ) =0 . 

DO 21 K=1 , M 

21 AX(IT)=AX(IT)+BIIT,K)*A(K,J) 

36 01 J)=-AXI IT) 

I NDE X = 0 

GO TO 8 
END 
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SUBROUTINE ELEM 



IMPLICIT REAL*8 <A-H,0-Z> 

COMMON A (70, 150) , B<70,70) ,C< 150 ) t D< 1501 ,H<70) ,AX< 701 
COMMON DB ( 70 ) » DUAL ( 70 >,R<150),CE(150),E,W,Z,DI,FF,EP 
COMMON CNAMEI 150) 

COMMON IB( 70) ,M,N, JS, I R , I NOE X , I TR , Ml , N1 
I R=0 
J= JS 

CALCULATING THE REVISED SIMPLEX TABLEAU PIVOT COLUMN 
WITH ROUNO-OFF ERROR CONTROL 

DO 22 1=1, M 
DM=0. 

DP = 0 « 

AX ( I )=0. 

DO 21 K= 1 , M 
T=B ( I , K ) *A ( K , J ) 

IFIT.GE.O. ) GO TO 40 
DM=DM+T 
GO TO 21 
40 DP=DP+ T 

21 CONTINUE 

AX ( I )=DP+DM 
I F ( DP . EQ. 0 • ) GO TO 22 
IF<DABS<AX< I ) /DP) .GT.EP) GO TO 22 
AX ( I ) =0 • 

22 CONTINUE 

FINDING THE PIVOT ROW 

JR=0 
EE = E 

31 DO 25 1=1, M 

I F < INDEX. EQ.O) GO TO 5 
IF ( I B ( I ) ) 5,25,25 
5 IF(AX< I )-E) 25,25,23 

23 IF( IR.EO.O) GO TO 24 

IF (DABS (HI I ) -AX ( I ) *X ) -EE ) 70, 70, 2 

2 I F < H < I)-AX< I >*X >24,70,25 

24 X=H ( I ) / AX ( I ) 

3 IR= I 

EE=E / AX ( IR) 

GO TO 25 

70 IF t I 8 < I >> 74,25,73 

73 IF < I B < IR) )25,25,4 

74 IF( IB< IR) >4,25,3 

4 I F ( A X ( I R ) -AX ( I ) ) 25 , 25 , 3 

25 CONTINUE 
RETURN 
END 
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c 

c 



SUBROUTINE PIVOT 



C 

C 

C 



IMPLICIT REAL*8 (A-H,0-Z) 

COMMON A<70, 150) ,B{ 70,70) ,C< 150 ) ,0< 150) ,H(70) ,AX< 70) 
COMMON DB ( 70 ) , DUAL ( 70 ) , R ( 1 50 ) , C F ( 1 50 ) , F , W , Z , 0 1 , FF f EP 
COMMON CNAME ( 150) 

COMMON IB(70),M,N,JS,IR, INDEX, ITR,M1,N1 
DIMENSION P ( 200 ) , BX ( 200 ) 

WRITE (6,214) IR, AX(IR) 

214 FORMAT ( 1H, 25X, 'IR=', 13, 3X, •AX=», F12.4) 

PERFORMANCE OF THE PIVOT 



43 IF ( I B ( IR ) ) 1,1,2 

2 K=IB< IR) 

R ( K ) =0 • 

GO TO 3 

1 P(JS)=1. 

3 I B ( I R ) = JS 

DI=D I *DABS( AX ( IR) ) 

IF(DI,GE. 9999998.) GO TO 40 
DI=IDINT(DI+0.5) 

40 CONTINUE 
E=0. 5/DI 
DO 27 1=1, M 

27 P( I )=-AX< I ) /AX ( IR) 

P( I R ) =1 ♦ / AX ( I R ) - 1 . 

P ( M+l ) =-D ( J S ) /AX ( IR) 

DO 28 J= 1 , M 
BX ( J ) =B ( I R , J ) 

28 CONTINUE 
HX=H ( IR ) 

DO 29 J= 1 , M 
DO 29 1=1, M 
T=B ( I , J) 

B( I , J)=B( I , J) +P( I ) *BX ( J ) 

IF(T .EO.O. ) GO TO 29 
IF(DABS(B(I ,J)/T).GT.EP) GO TO 29 
B( I , J) =0. 

29 CONTINUE 

DO 30 1=1, M 
T=H< I ) 

H< I) =H< I ) +P( I ) *HX 
I F ( T .EO.O. ) GO TO 30 
IF(DABS(H( I)/T) .GT.EP) GO TO 30 
H ( I ) =0 . 

30 CONTINUE 
W=W-P(Mfl)*HX 
DO 31 J=1,M 
T=DB ( J ) 

DB( J)=DB(J)+P(M+1) *BX ( J ) 
IF(T.EO.O. ) GO TO 31 
I F(OABS(DB( J) /T) .GT.EP ) GO TO 31 
DB ( J ) =0 • 

31 CONTINUE 



UPDATING THE DB, B, AND H VALUES 

IFCDI. GE. 999998. ) GO TO 81 
DO 92 J=1,M 
DF=DI*DABS(DQ(J))+0.5 
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I 



non 



IF(DF. GT. 999998. ) GO TO 61 
DF= I D I NT ( OF ) /0 I 
DB(J)=OSIGN(DF,DB(J ) » 

61 00 92 1=1, M 

0F=D I *D ABS ( B ( I , J ) ) +0. 5 
IF(DF. GT. 999998. ) GO TO 92 
DF= I D I NT ( OF ) /OI 
6( I »J)=DSIGN(OF,B( I,J ) ) 

92 CONTINUE 

00 93 1 = 1, M 

DF=D I *DABS ( H ( I ) ) +0 . 5 

IF(OF. GT. 999998. ) GO TO 93 

CF=IDINT(DF)/OI 

H ( I )=DSIGN< OF ,H ( I ) ) 

93 CONTINUE 
81 CONTINUE 

DF=W*DI 

IF(DF. GE. 9999998. ) GO TO 87 
W=IDINT(DI*W+0.5)/DI 
87 CONTINUE 

CALCULATION OF 0 ( J ) • S FOR NEXT ITERATION 

DO 39 J = 1 » N 
DM=0. 

DP=0 • 

0 ( J ) =0 . 

DO 33 K=1,M 
T=DB(K)*A(K,J) 

IFIT.GE.O.) GO TO 91 
0M=DM+T 
GO TO 33 
91 DP=DP+T 
33 CONTINUE 
0 ( J ) =DP + DM 

IFIDP.EO.O. ) GO TO 39 
IF(DABS(D(J)/OP).LE.EP) GO TO 35 
IF(DABS(D(J) )-E)35,35,39 
35 D( J) =0. 

39 CONTINUE 
0( JS )=0. 

I F ( W . GT .0.001) GO TO 99 
W=0 

99 I TR= I TR+1 

WRITE (6,210) ITR, W, Z, 01 
210 FORMAT ( 1H , 5X , 'ITR=», 15, 5X, • W= • , F18.8, 

1 5 X , *Z=', F12.9, 5X , *01 = ', E12.9) 

20 RETURN 

21 STOP 
FNO 
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